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We  present  a  new  multidomain  spectral  collocation  method  that  uses  staggered  grids  for 
the  solution  of  compressible  flow  problems.  The  solution  unknowns  are  defined  at  the  nodes 
of  a  Gauss  quadrature  rule.  The  fluxes  are  evaluated  at  the  nodes  of  a  Gauss- Lobatto  rule. 
The  method  is  conservative,  free-stream  preserving  and  exponentially  accurate.  A  significant 
advantage  of  the  method  is  that  subdomain  corners  are  not  included  in  the  approximation, 
making  solutions  in  complex  geometries  easier  to  compute. 
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1  Introduction 

Standard  Chebyshev  spectral  methods  applied  to  compressible  flow  problems 
have  some  severe  restrictions  [7].  The  computational  domain  must  be  simple  enough  to 
map  onto  a  square,  in  two  space  dimensions,  or  a  cube  in  three.  To  increase  spatial 
resolution  the  polynomial  approximation  order  must  be  increased.  For  high  orders,  the 
derivative  approximations  must  be  performed  with  Fast  Fourier  Transform  methods  to  be 
efficient.  If  matrix  multiplication  is  used  instead,  the  work  grows  too  rapidly  with  the 
number  of  degrees  of  freedom  to  be  practical.  Finally,  the  time  step  restrictions  are 
severe  since  the  time  step  decreases  asymptotically  as  the  square  of  the  order  of  the 
approximating  polynomials. 

The  basic  premise  of  a  multidomain  method  is  that  these  restrictions  can  be 
reduced  by  subdividing  the  computational  domain  into  multiple  zones,  called 
subdomains,  on  which  the  spectral  approximation  is  applied.  As  a  result,  the  method  can 
be  used  on  more  complex  geometries.  The  use  of  lower  order  approximating 
polynomials  in  each  subdomain  means  that  matrix  multiplication  can  be  both  efficient 
and  accurate,  and  the  time  step  restrictions  need  not  be  as  severe.  A  discussion  of  the 
advantages  of  multidomain  methods  over  the  single  domain  method  was  presented  in  [24] 
and  has  been  updated  in  [19]. 

Less  than  a  decade  after  they  were  first  introduced[24],  the  bulk  of  the  spectral 
multidomain  methods  for  that  have  been  proposed  compressible  flows  or  similar 
hyperbolic  problems  still  define  the  solution  unknowns  at  the  nodes  of  the  Gauss-Lobatto 
quadrature,  just  as  in  a  single  domain  method.  Examples  include  [33],  [25],  [29]  and  [3] 
for  general  hyperbolic  problems,  and  [26]  for  the  Euler  gas-dynamics  equations. 
Methods  for  the  advective  terms  of  the  compressible  Navier-stokes  equations  were 
presented  in  [18]  and  [19],  An  interesting  method  for  coupled  acoustic  and  elastic  wave 
interactions  was  proposed  in  [1]. 


1 


The  main  differences  between  the  Lobatto  grid  methods  are  whether  the  equations 
are  written  in  conservative  or  non-conservative  form,  and  the  manner  in  which  the 
interfaces  are  treated.  The  conservative  form  of  the  equations  was  used,  for  example,  in 
the  methods  presented  by  [18]  and[4].  Non-conservative  forms  were  considered  in  [33], 
[26]  and  [1].  We  will  show  below  that  the  use  of  the  conservative  form  of  the  equations 
does  not  guarantee  that  the  method  is  globally  conservative,  since  the  interface  treatment 
may  lead  to  loss  of  conservation. 

Two  approaches  have  been  used  at  subdomain  interfaces  to  ensure  that  waves 
propagate  properly  through  them.  The  two  methods  were  contrasted  in  [25]  and  are 
considered  in  more  detail  in  [3].  At  least  two  values  of  the  normal  derivative  are 
available  an  interface  point,  depending  on  the  number  of  subdomains  that  have  that  point 
in  common.  One  interface  method  integrates  a  differential  compatibility  equation  for  the 
points  along  the  interface  [24],  [33],  [18],  [3].  Derivatives  are  chosen  from  appropriate 
subdomains  so  that  wave  components  are  “upwinded”.  The  other  approach  uses  a 
correction  procedure  ([25],  [26],  [1],  [4]).  To  implement  the  correction  method,  the 
interior  point  approximation  is  integrated  everywhere,  including  at  the  boundaries.  As  a 
result,  multiple  solution  values  are  available  at  each  interface  point.  A  characteristic 
combination  of  these  solutions  is  then  made  to  correct  the  solution  for  the  propagation  of 
waves  across  the  interface. 

Each  interface  treatment  has  its  advantages  and  disadvantages.  Integrating  the 
differential  equation  means  that  the  solution  can  be  approximated  to  any  order  accuracy 
in  time,  depending  on  the  choice  of  the  time  integration  scheme.  Implicit  time  integration 
schemes  can  also  be  used[3].  The  serious  disadvantage  is  that  the  tangential  space 
derivatives  must  be  continuous  across  subdomain  interfaces  in  more  than  one  space 
dimension.  This  requirement  severely  restricts  the  types  of  geometries  on  which  solutions 
can  be  computed,  since  it  means  that  the  Jacobians  of  the  transformations  of  the 
mappings  between  the  subdomains  and  the  unit  square  must  be  continuous  across 
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subdomain  interfaces.  The  correction  scheme,  on  the  other  hand,  does  not  require 
smoothness  of  the  grids,  since  only  the  solution  values  are  required  at  the  interfaces. 
However,  the  temporal  accuracy  of  the  correction  method  is  limited  to  first  order.  (C.f. 
[7],  page  245.) 

A  disadvantage  shared  by  the  two  interface  treatments  is  their  complexity.  Either 
method  is  simple  to  apply  in  one  space  dimension.  In  two  space  dimensions  a  choice 
must  be  made  at  corners  to  determine  from  which  subdomains  the  solution  bi¬ 
characteristics  must  be  computed.  Special  algorithms  can  be  developed  for  the 
approximations  at  the  corners  of  subdomains  [26],  but  if  more  than  four  subdomains  meet 
at  a  single  point,  the  choice  can  be  even  more  complex. 

A  very  different  multidomain  approximation  is  based  on  the  Chebyshev  cell 
averaged  grid  originally  proposed  by  Cai  et  al.  [5].  In  this  method,  “cell”  averaged 
quantities  are  defined  on  the  Gauss-Chebyshev  grid,  while  fluxes  are  defined  at  the  more 
usual  Lobatto  points  [35],  [15],  [16].  The  cell  averaged  method  avoids  many  of  the 
disadvantages  of  the  methods  just  described.  It  is  fully  conservative,  and  it  can  be 
approximated  to  any  temporal  order  of  accuracy.  It  is  also  geometrically  flexible  because 
it  does  not  require  continuity  of  the  transformations  across  interfaces.  In  more  than  one 
space  dimension,  the  method  does  require  special  attention  at  the  corners  of  subdomains. 
Currently,  a  simple  average  of  the  multiple  solutions  is  computed  and  broadcast  to  all 
contributing  subdomains  [16]. 

In  this  paper,  we  present  a  new  multidomain  spectral  collocation  method  for  the 
solution  of  compressible  flow  problems.  The  new  method  is  based  on  a  staggered  grid, 
analogous  to  fully  staggered  grids  often  used  with  finite  difference  methods.  The 
solutions  are  defined  at  the  nodes  of  a  Gauss  quadrature  rule,  and  the  fluxes  are  evaluated 
at  the  nodes  of  a  Gauss-Lobatto  rule.  Staggered-grid  spectral  approximations  were  first 
proposed  for  the  solution  of  the  incompressible  Navier-Stokes  equations.  (C.f.  [7],  page 
234.)  Our  grid  will  be  identical  to  the  fully  staggered  grid  of  Bemardi  and  Maday  [2]. 
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The  staggered  grid  multidomain  method  for  compressible  flow  problems  has  all 
the  desirable  features  found  in  the  methods  discussed  above.  First,  like  the  cell  averaged 
method  it  is  conservative.  Thus,  it  should  be  possible  to  apply  shock  capturing 
techniques  to  the  approximation.  Subdomains  can  be  defined  independently  of  their 
neighbors,  so  the  method  is  geometrically  flexible.  The  interface  condition  can  be 
computed  to  the  same  temporal  accuracy  as  the  interiors.  However,  in  multiple  space 
dimensions  the  method  does  not  include  (the  Gauss  rules  being  open)  the  corners  of 
subdomains.  Thus  the  coding  of  the  method  does  not  require  special  cases  at  corners,  and 
any  number  of  subdomains  can  meet  at  a  point  without  difficulty. 

The  paper  is  divided  as  follows.  The  algorithm  is  presented  in  the  next  section  for 
problems  in  one  space  dimension,  along  with  definitions  of  the  notation  used  throughout 
the  paper.  We  show  that  the  staggered  grid  method  is  conservative,  while  methods  that 
upwind  derivatives  are  not.  A  scalar  problem  and  a  linear  system  will  be  used  as 
examples  to  show  that  the  method  is  exponentially  convergent  for  smooth  problems. 
Though  we  will  be  concerned  in  this  paper  primarily  with  steady  problems,  an  example  is 
included  to  show  that  high  order  temporal  accuracy  can  also  be  obtained.  In  Section  3, 
we  describe  the  algorithm  in  two  space  dimensions.  We  show  that  the  method  remains 
conservative  and  is  also  free-stream  preserving.  Section  4  provides  three  examples  of  the 
use  of  the  method  for  two-dimensional  problems.  The  first  problem  is  that  of  a  point 
source  flow,  for  which  there  is  an  exact  solution.  We  show  that  exponential  accuracy  is 
obtained  for  this  problem.  The  second  problem  is  a  subsonic  flow  over  a  circular  bump 
in  a  channel,  and  we  show  that  the  entropy  errors  decay  exponentially  fast.  Finally,  we 
solve  a  transonic  flow  in  an  axisymmetric  converging-diverging  nozzle  and  compare  the 
results  to  experimental  data.  Concluding  remarks  are  then  made  in  the  Section  5. 


4 


2  The  Staggered  Grid  Approximation  in  One  Space 
Dimension. 

2.1  Notation 

The  staggered  grid  approximation  uses  two  grids  to  compute  the  solution  values 
and  advective  fluxes.  Unlike  the  common  Chebyshev  approximation  [7],  which  uses 
only  the  nodes  of  the  Gauss-Lobatto  quadrature  as  collocation  points,  the  new  method 
uses  both  the  Gauss  and  the  Gauss-Lobatto  points.  We  denote  the  points  on  the  two  grids 
by  the  Lobatto  points,  Xj,  and  the  Gauss  points, 

.  .  (1) 
1\  \2N  +  2  JJ 

In  (1),  we  have  mapped  the  usual  collocation  points  defined  on  [-1,1]  to  the  more 
convenient  unit  interval.  The  overbar  and  half  point  notations  for  the  Gauss  points  are 
used  only  for  their  value  as  an  analogy  to  staggered  grid  finite  difference  methods.  It 
must  be  understood  that  the  Gauss  points  do  not  lie  halfway  between  the  Lobatto  points 
[7]. 

Two  polynomial  approximations  are  defined,  one  for  each  grid.  Let  the  space  of 
polynomials  of  degree  less  than  or  equal  to  N  be  denoted  P^.  Let  £j(i^)€P^  be  the 

Lagrange  interpolating  polynomial 

A-«)=fl 

i=0 

defined  on  the  Lobatto  grid.  On  the  Gauss  grid,  we  define  hj^y2  ^  be  the 

polynomial 

(2b) 

1=9  -^i+i/2 
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Finally,  let  Qj  be  a  grid  point  value  on  the  Lobatto  grid  and  be  a  value  defined  on 
the  Gauss  grid.  Then  we  write  the  polynomials  that  interpolate  these  values  as 


N 


Q(X)  =  J^Q/j(X) 

j=0 

(3a) 

(3b) 

j=0 


2.2  The  One  -Dimensional  Staggered  Grid  Approximation  for  Scalar  Equations 

To  motivate  the  staggered  grid  approximation,  we  consider  the  approximation  of 
scalar  problems  of  the  form 

u, +/^(i<)  =  0  df ldu>Q,x&[a,b\,t>Q 
<  u(x,Q)  ^  Uf^(x)  (4) 

u(a,t)  =  g(t) 

The  interval  [a,b]  is  subdivided  into  multiple,  non-overlapping  subdomains,  =[cik,h]^ 
k  =  which  are  ordered  left  to  right.  A  simple  linear  transformation  can  be  made 

to  the  unit  interval,  so  that  on  each  subdomain  we  solve  the  problem 

u,+—f,W  =  0  Xe[0,l],t>0  (5) 

x^ 

On  each  subdomain  is  placed  the  staggered  grid  defined  by  (1).  For  convenience, 
we  will  assume  that  the  same  number  of  points  is  used  in  each  subdomain,  but  this  is  not 
required  by  the  method.  We  then  let  defined  by  (3b),  approximate  the 

exact  solution,  u  on  Similarly,  the  flux  is  approximated  by  the  polynomial 
F'^iX)  €  P^,  defined  by  (3a).  Substitution  of  these  approximations  into  (5)  gives 

=  k  =  l,2,...,K  (6) 

x^  oX 
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To  obtain  the  equations  that  define  the  solution  unknowns  at  the  Gauss  points,  we  require 
that  the  residual,  R,  be  zero  at  the  Gauss  points  of  the  subdomain.  This  leads  to  the 
collocation  approximation 


(7) 


dt  dX 

The  spatial  derivative  operation  in  (7)  can  be  evaluated  as  the  multiplication  of  the 
vector  of  flux  values  by  a  derivative  matrix,  D.  From  (3a),  we  see  that 


ax 


n=0 


n=0 


(8) 


so  we  write 


dX 


=  (DFM  =^d-Fl 

V  /;+l/2  1"  ” 


;+!/2 


n=0 


(9) 


Thus,  (7)  can  be  written  in  vector  form  as 

— +  DF'^=0  k  =  l,2,...,K  (10) 

dt 

where  U*  =  [t/f„  C/,‘,, ...  C/J.,„]''.F‘  =  [fJ  Ff  ...  Fjf . 

To  compute  the  flux  values  on  the  Lobatto  grid,  we  use  the  following 
reconstruction  procedure.  We  first  evaluate  the  interpolant  U’‘{X)  e  P^_i  at  the  Lobatto 

points  by  multiplying  the  vector  of  solution  values  by  an  interpolation  matrix,  I. 

£7(X,)  =  =  (11) 

n=0  «=0 

The  family  of  characteristics  of  (5)  runs  left  to  right.  Thus,  we  expect  the  use  of  the 
solution  extrapolated  to  the  left  subdomain  boundary  to  lead  to  an  unstable  procedure.  To 
provide  the  proper  characteristic  domain  of  dependence,  we  use  the  boundary  condition 
to  define  the  7  =  0  value  on  the  furthest  left  subdomain.  At  subdomain  interfaces,  where 
two  values  (0)  are  available,  we  choose  the  value  computed  from  the  left  side 
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of  the  interface.  The  result  is  an  upwind  evaluated  approximation  at  both  the  left 
boundary  and  the  interfaces.  The  fluxes,  Fj,  are  then  computed  from  the  solution,  values 
on  the  Lobatto  grid. 

The  method  imposes  the  boundary  conditions,  weakly,  through  the  definition  of 
the  flux,  since  the  discrete  solution  values  are  not  used  directly  at  the  boundary  or 
interfaces.  To  see  this,  consider  the  single  domain  approximation  of  (4)  for/=  u.  Then 
we  can  write  the  flux  F(Z)  =  f/(X)  e  in  terms  of  the  interpolant  TJ{X)  and  the 

boundary  condition  as 


V{X)^U{X)  +  [g-U{a)\l,{X),  (12) 


so  that  the  polynomial  U(X}  satisfies 


U{X.)  = 


\U{Xj)  7  =  1,2, 

\g  7  =  0 


Then  (7)  can  be  written  as 


N 


(13) 


(f/'(x„)-«]C(X-.«)  >=0.1 . Af-i-  (14) 

Cl-t 

Thus,  the  boundary  condition  is  imposed  indirectly  at  each  collocation  point  through  the 
penalty  term  on  the  right. 

Equation  (10)  is  a  system  of  ordinary  differential  equations  that  must  be 
integrated  in  time  to  get  the  approximate  solution  values  at  the  Gauss  points.  In  principle, 
any  common  integration  procedure  can  be  used.  We  have  chosen  to  use  low  storage 
Runge-Kutta  methods  that  require  only  2-N  storage  locations.  For  the  computation  of 
steady-state  problems,  for  which  the  time  discretization  is  only  an  iterative  procedure,  we 
have  used  a  mid-point  rule. 


j  j  k,n+\/2  _  _ 

^7+1/2  ~  ^7+1/2 


Ar  1  dF 


<k,n 


2  dX 


;+l/2 


^  j+in 


1  dF 


'k.n+l/2 


dx 


7  =  0,1,...,7V-1 


7  =  0,l,...,iV-l 


j+l/2 


(15) 
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This  method  appears  to  have  a  good  balance  between  the  time  step  that  can  be  used  and 
temporal  damping  introduced  by  the  scheme.  With  additional  knowledge  of  the 
eigenvalue  structure  of  the  differentiation  matrices,  other  choices  might  include  schemes 
optimized  for  rapid  convergence  to  steady-state,  such  as  those  discussed  in  [11],  [10]  and 
[9]. 

For  time  dependent  problems,  we  have  used  the  third  order  2-N  storage  method  of 
Williamson  [39],  and  the  more  recent  fourth  order  scheme  of  Carpenter  and  Kennedy  [8]. 
We  note  that  it  should  also  be  advantageous  to  use  the  new  low  storage  Runge-Kutta 
methods  derived  by  Hu  et  al.  [21],  which  are  optimized  to  minimize  the  phase  and 
dissipation  error  introduced  by  the  temporal  approximation. 

To  summarize  the  staggered  grid  procedure,  we  present  the  following  algorithm 
for  the  scalar  problem  described  above: 

Algorithm  I.  (Staggered  Grid,  Scalar,  ID) 

1.  Interpolate  U  =  [Ui/2’U3/2>-"’UN-i/2r  Lobatto  points: 

Compute  the  matrix- vector  product  U*"  =  defined  by  (1 1)  for  each  subdomain 

2.  Compute  the  flux  values  at  internal  points  on  Lobatto  Grid: 

F]=f{U])  y  =  i,2,...,yv,k  =  i,2,...,K 

3.  Apply  boundary  and  interface  conditions: 

^0  =  f(.8) 

Fo=f{K')  k^2,3,...,K 

4 .  Differentiate  Flux  and  evaluate  on  Gauss  grid  by  subdomain: 

Compute  the  matrix- vector  product  k  =  l,2,...,A^by  eq.  (9). 

5.  Update  the  solution  by  subdomain: 

Integrate  (10)  by  the  chosen  ODE  solver,  repeating  Steps  1-4,  as  necessary. 

6.  Repeat  1-5  until  done. 
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We  note  that  the  method  requires  two  matrix- vector  products  per  Runge-Kutta 
stage.  This  is  twice  the  work  of  a  Lobatto  grid  method,  or  of  the  cell-averaged  method. 
Thus,  there  is  no  speed  advantage  for  the  method  in  one  space  dimension. 

A  desirable  feature  of  the  staggered  grid  approximation,  (7),  is  that  it  is 
conservative.  To  show  conservation,  we  define  the  quadrature 


N-l 


0  ;=o 

1 

2=lh,,m(X}dX 


(16) 


i+1/2 


For  each  j,  we  multiply  (7)  by  The  sum  over  all  points  and  all  subdomains  is 


K  N-l 

IE 

i!:=l  ;=0 


W 


i+l/2 


dU; 


dt 


j+Ml 


=  0. 


(17) 


Now,  U  (X),F  (X)  e  and  is  a  polynomial  of  degree  zero,  so  we  can  replace  the 
sum  over  j  by  integrals  to  get 


I 

SI 

k=l 


dt 


V  =  0. 


(18) 


Upon  integrating  the  flux  derivatives,  the  interface  contributions  cancel,  and 


-f  X  \v\^)4dx\  =  F‘(0) -  F^(l)  (19) 

dt  J 

In  contrast,  a  multidomain  method  defined  on  the  Lobatto  points,  where  the 
upwind  value  of  the  flux  derivative  is  used  at  the  interface  (e.g.,  [18])  is  not  truly 
conservative.  To  show  this,  it  is  sufficient  to  consider  two  subdomains,  and  for 
which  =  1.  Using  the  upwind  derivatives  at  the  interface,  one  integrates  the  following 
system  of  equations 
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Vo  =8 
dU^ 

-^  =  -F'V  j  = 
dt 

dUf 
dt 


-F;"  ;  =  l,2,...,iV-l 


dUl^  _  dU, 


(20) 


0  _ 


dt 


dt 


-f: 


./i 


p.. 

We  then  multiply  each  equation  by  its  associated  Clenshaw-Curtis  weight  [7]  and  sum  to 
get 

^.dUf  ^  ^dUf  , 

y — ^wf+y — ^wf  = 

Po  dt  ^  dt  ^ 


w. 


'"o'" + -  -f.'"] 

Ut  J  j^Q 


(21) 


=  FHO)  -  F\l)  -  Wo"[n,"  -  Fo'"]  +  Wo" 


^  dt 


Eq.  (21)  shows  that  there  is  a  contribution  at  the  interface  proportional  to  the  jump  in  the 
derivative  across  the  interface.  For  smooth  enough  functions,  this  is  not  expected  to  be  a 
problem,  since  the  difference  between  the  derivatives  should  go  to  zero  exponentially 
fast,  while  the  coefficient  decays  as  0{llN^). 

As  an  example  of  solutions  computed  using  Algorithm  I,  we  compute  a  steady 
solution  of  the  equation  u^  +  u^=  f  x  e  [0,2],t  >  0 .  Scalar  examples  of  one  dimensional 

time  dependent  problems  can  be  found  in  [27].  Comparisons  to  a  variety  of  finite 
difference  methods  can  be  found  in  [37].  The  initial  and  boundary  conditions  were 
chosen  so  that  the  exact  steady  solution  is  u{x)  =  tanh((x-1.5)/  2).  Fig.  1  shows  the 
solution  computed  using  three  subdomains  and  eight  points  per  subdomain. 
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Fig.  1.  Steady  solution  of  a  scalar  wave  equation 

Convergence  of  the  error  is  exponential  for  this  problem.  Fig.  2  shows  the  error 
plotted  as  a  function  of  the  Gauss  polynomial  order  for  the  subdivision  shown  in  Fig.l. 
For  comparison,  we  have  also  plotted  the  error  of  the  single  grid  multidomain  method 
(20).  We  see  that  the  staggered  grid  approximation  is  at  least  as  accurate  as  the  non- 
staggered  grid  approximation,  and  sometimes  more  accurate  by  a  factor  of  four. 
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Polynomial  Order 

Fig.  2.  Convergence  of  the  error  for  the  staggered  grid  multidomain  method  compared  to 
a  non-staggered  multidomain  approximation. 

2.3  The  One  -Dimensional  Staggered  Grid  Approximation  for  Systems 

Algorithm  I  can  be  easily  extended  to  systems  of  hyperbolic  equations  of  the  form 

|Q.  +  F,(Q)  =  0  xe[aMt>0 
|q(x,0)  =  Qo(x) 

where  Q  and  F  are  m-vectors.  We  assume  that  the  system  is  hyperbolic,  that  is,  the 
Jacobian  matrix  A  =  dF  /dQ  =  ZAZ“\  where  A  is  a  real  diagonal  matrix.  We  further 
assume,  as  is  the  case  for  the  Euler  gas-dynamics  equations,  that  the  flux  can  be  written 
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as  F  -AQ.  To  complete  (22),  we  assume  that  appropriate  dissipative  boundary 
conditions  are  applied. 

The  approximation  of  the  system  follows  that  of  the  scalar  equation,  except  for 
the  treatment  of  boundary  and  interface  conditions.  At  an  interface  between  subdomains 
k-l  and  k,  there  are  two  vector  values  of  the  interpolated  solution  available,  and  Q*. 

The  computed  flux  must  use  these  two  values  to  allow  waves  to  propagate  through  the 
interfaces.  For  constant  coefficient  linear  problems,  we  can  write 

F  =  AQ  =  ZAZ 'Q  =  ZA^Z 'Q  +  ZA  Z 'Q  (23) 

where  A“  =  A±jA|.  The  first  term  represents  waves  moving  left  to  right,  and  the  second 
represents  waves  moving  right  to  left.  An  upwind  approximation  chooses  for  the 
right  going  components,  and  Q*  for  the  left  going  components  to  give 

=  K  =  !F{Qi7' ,  Q9  =  ZA"Z-'Q^^-'  +  ZA'Z-'Q^ .  (24) 

Characteristic  decompositions  for  nonlinear  flux  vectors  have  been  addressed 
extensively  in  the  finite  difference  community  (Ref.  [20]).  We  have  considered  flux 
vector  splitting  and  flux  difference  splitting. 

The  resolution  of  the  jump  at  the  subdomain  interfaces  can  be  easily  viewed  using 
flux  vector  splitting.  The  flux  is  decomposed  into  a  right  going  and  a  left  going  flux, 
F(Q)  =  F'^(Q)  +  F  (Q).  The  splitting  is  done  so  that  the  Jacobian  matrix  of  F+  has  only 
positive  eigenvalues  and  the  Jacobian  of  F‘  has  only  negative  eigenvalues.  Examples  are 
the  Van  Leer  [38]  splitting  and  the  more  recent  splitting  of  Liou  and  Steffen  [30].  Using 
flux  vector  splitting,  the  positive  flux  is  evaluated  using  the  solution  from  the  left,  the 
negative  flux  is  computed  using  the  value  from  the  right 

=  F"(Q^^-')  +  F-(Q^o)-  (25) 

Van  Leer’s  Flux  vector  splitting  was  used  in  [18]  to  compute  the  flux  derivatives  at 
interfaces  for  the  advective  part  of  the  Navier-Stokes  equations. 
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As  an  interface  treatment,  flux  vector  splitting  has  the  desirable  feature  that  the 
positive  and  negative  fluxes  can  be  computed  within  a  subdomain  without  regard  to  the 
neighbors.  The  final  flux  computation,  (25),  requires  only  a  simple  sum  of  the  boundary 
fluxes.  However,  we  found  that  the  Van  Leer  splitting  applied  to  the  staggered  grid 
scheme  was  unstable  for  some  long  time  integrations. 

As  an  alternative  method  to  compute  the  interface  flux,  we  have  chosen  to  use  an 
approximate  Riemann  solver.  This  approach  was  also  taken  by  Giannakouros  and 
Karniadakis[15].  Several  solver  choices  are  possible,  but  we  have  used  Roe’s  [34]  solver 
with  the  entropy  fix.  Formally,  given  the  two  states  and  Qq,  we  write 

y{Qjr',Q;)=i(F(Q‘„-')+F(Q5))-^i«|Ar'(Q:  -Q»‘)  <26) 

where  R  is  the  matrix  of  the  right  eigenvectors  of  the  Jacobian  of  F,  computed  using  the 
Rop-average  of  and  Qq.  The  eq.  (26)  is  modified  to  correct  the  entropy  across  sonic 

points  [20]. 

Boundaries  can  be  considered  as  interfaces  between  the  computed  solution  and 
the  solution  assumed  to  exist  outside  the  computational  region,  if  fully  known.  Thus,  we 
can  compute  the  boundary  flux  by 


Fl=!F{Q(a,t),Ql) 


(27a) 


on  the  left,  and 

=  ir(Q^,Q(^,0)  (27b) 

on  the  right,  where  Q(a,  t)  and  Q{b,  t)  represent  the  exterior  solution  at  the  boundaries. 
Other  ways  to  compute  the  boundary  flux  when  the  full  exterior  solution  is  not  known 
will  be  described  in  regard  to  specific  problems  in  Section  4. 

In  summary,  for  systems  of  equations,  we  have  the  following  algorithm: 
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Algorithm n.  (Staggered  Grid,  System,  ID) 


1 . Interpolate  the  Gauss-point  solution  values  to  the  Lobatto 
points : 

Compute  the  matrix- vector  product  by  eq.  (11)  for  each  subdomain. 

2.  Compute  the  interior  point  fluxes: 

F;=F(q;)  j  = 

3.  Apply  the  interface  conditions: 

F*^-'=F^=y-(Q^-^Q^)  k  =  2,...,K 
3.  Apply  boundary  conditions  at  left  and  right 


4.  Compute  spatial  derivatives  at  Gauss  points: 

Compute  the  matrix-vector  product  by  eq.  (9). 

5.  Update  the  solution  at  the  Gauss  points 


+  F'1u2=0  j  =  0,l,...N-l,k  =  l,2,...,K, 


dt  - 

repeating  Steps  1-4  for  each  Runge-Kutta  stage. 

6.  Repeat  Steps  1-5  until  done 


As  an  example  of  the  application  of  Algorithm  II,  we  solve  the  problem 
Q,  +  F^=0  jc6[-l,4],r>0 


Q  = 


(28) 


Further  examples  that  model  unsteady  acoustic  propagation,  including  acoustic 
propagation  in  a  quasi-one-dimensional  nozzle,  can  be  found  in  [27].  The  initial 
condition  for  (28)  was  chosen  to  be  a  Gaussian  pulse  with  the  peak  at  x  =  1 


Q(x,0) 


0 


(29) 


We  specify  boundary  conditions  so  that  the  waves  pass  through  the  boundaries  without 
reflection. 
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(30) 


m(4,0  +  v(4,0  = 

The  solution  was  computed  using  four  subdomains  of  equal  length  and  with  the 
same  number  of  points  within  each  subdomain.  For  the  time  discretization,  we  have  used 
both  the  Williamson  [39]  order  and  the  Carpenter  and  Kennedy  [8]  4'^  order  schemes. 
Fig.  3  shows  the  solution  at  time  t  -  0.75.  For  a  small  enough  time  step,  spectral  decay 
of  the  error  is  observed,  as  shown  in  Fig.  4.  To  study  the  temporal  accuracy,  we 
computed  the  solution  on  the  finest  grid,  N  =  25,  so  that  the  spatial  accuracy  was  close  to 
rounding  error.  A  plot  of  the  errors  as  a  function  of  At  for  the  third  and  fourth  order 
methods  is  shown  in  Fig.  5.  A  least  squares  fit  to  the  errors  indicates  a  slope  of  2.995  for 
the  third  order  and  3.998  for  the  fourth  order  method,  so  the  expected  high  order  temporal 
accuracy  is  obtained.  We  also  see  that  over  the  At  range  shown,  the  fourth  order  method 
is  about  two  orders  of  magnitude  more  accurate  than  the  third.  This  is  consistent  with  the 
observations  of  [8]. 


Fig.  3.  Solution  of  the  system  (27)  at  t  =  0.75  using  four  subdomains 
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Polynomial  Order 


Fig.  4.  Decay  of  the  spatial  error  of  the  system  (28). 


3  The  Two-dimensional  Approximation. 

We  now  describe  the  approximation  of  the  Euler  equations  of  gas-dynamics  in 
conservative  form, 


dQ  dF  dG  - 
— -I-  —  +  —  =  0, 
dt  dy 


(31a) 


where  Q  is  the  vector  of  solution  unknowns  and  F(Q)  and  G(Q)  are  the  advective  flux 
vectors 


pu 

pv 

pu 

F  = 

p  -t-  plp 

G  = 

puv 

.  2 

pv 

puv 

p  +  pv 

u(pe  +  p) 

_v(pe  +  p)_ 

(31b) 


We  assume  y  =  1 .4  and  that  pe  =  pf  (y-l)  +  (u^  +v^)/2.  For  axisymmetric  problems, 


such  as  the  transonic  flow  in  the  converging-diverging  nozzle  discussed  later,  we 
interpret  x  as  the  axial  coordinate  and  y  as  the  radial  coordinate.  We  then  add  to  the  right 
hand  side  of  (31a)  the  vector 


H 


pv 

puv 

0 

pv- 

v(pe  +  p) 


(32) 


3.1  Mapping  in  two  space  dimensions. 

In  two  space  dimensions,  we  subdivide  a  computational  domain,  O,  into 
quadrilateral  subdomains,  Q^,  k  =  Figure  6  shows  an  example  of  a  division  of  a 

region  into  four  subdomains.  We  make  three  assumptions  about  the  subdivision  in  this 
paper.  First,  we  allow  subdomains  to  intersect  only  at  a  point  or  along  an  entire  side. 
Second,  we  assume  that  the  approximation  is  conforming,  so  that  grid  lines  coincide 
across  subdomain  interfaces.  Finally,  we  assume  that  the  subdomain  boundaries  do  not 
move  in  time.  In  the  discussion  that  follows,  we  will  make  the  assumption  that  the  same 
polynomial  order  is  used  in  each  space  direction  and  for  each  subdomain.  In  practice,  the 
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number  grid  points  can  vary,  as  long  as  the  approximation  remains  conforming  at 
subdomain  interfaces. 


Fig  6.  Diagram  of  a  subdomain  decomposition  in  2D 


Subdomains  are  mapped  onto  the  unit  square  by  an  isoparametric  mapping.  Let 
the  vector  function  g(s),  0<  5  <1  define  a  parametric  curve.  Define  also  the  polynomial  of 
degree  N  that  interpolates  g  at  the  Gauss-Lobatto  points  to  be 

=  (33) 

j=0 

Four  such  polynomial  curves,  rin(5),  m  =  1,2,3,4,  counted  counter-clockwise,  bound  each 
subdomain.  We  map  each  subdomain  onto  the  unit  square  by  the  linear  blending  formula 

x"'  (X,  F)  =  (1  -  F)r,(Z)  +  Fr3(X)  -H  (1  -  X)r,  (F)  -f  xr^  (F) 
-Xj(1-X)(1-F)-X2X(1-F)-X3XF  ’ 

where  the  xfs  represent  the  locations  of  the  comers  of  the  subdomain,  counted  counter¬ 
clockwise. 

Under  the  mapping  ^  [0,l]x[0,l]  given  by  (34),  the  Euler  equations  (31) 
become 


dQ  ^ 

dt  J  dx^  dY 


=  0 


(35a) 


where 
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(35b) 


F  =  G  =  -y^F  +  x^G 

JiX,Y)  =  x^^y^-x^y^ 

Since  we  assume  here  that  the  subdomain  boundaries  do  not  move  in  time,  we  can  write 
(35a)  as 

^  +  ^  +  ^  =  0  (36) 

dt  dX  dY 

where  Q  =  /Q  and  the  fluxes  are  still  defined  by  (35b). 


3.2,  The  staggered  grid 

A  fully  staggered  grid  is  used  in  two  space  dimensions.  A  schematic  of  the  grid 
on  a  single  subdomain  is  shown  in  Fig.  7.  The  grid  is  the  same  as  the  staggered  grid 
proposed  by  Bemardi  and  Maday  [2]  for  the  solution  of  the  incompressible  Navier- Stokes 
equations.  In  what  follows,  we  will  ignore  superscripts  that  denote  which  subdomain  is 
being  considered,  unless  necessary. 
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Fig.  7.  Diagram  of  the  fully  staggered  grid  in  two  space  dimensions. 


Points  of  type  “a”  in  Fig.  7  represent  the  Gauss/Gauss  points 
(X;+i/2,F^+j/2),  1,7  =0,1,...,A-1.  The  grid  that  results  from  these  points  is  the  tensor 

product  of  the  one  dimensional  Gauss  grid  defined  in  (1).  We  approximate  the  solution 
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and  the  transformation  Jacobian  at  the  Gauss/Gauss  points,  and  denote  them  by 
Qi+i/2.;+i/2  and  =  ■^(^i+i/2,y+i/2’i^+i/2.;+i/2)-  Ffom  these,  we  compute  the  Gauss 

point  values  Q,+i/2,;+i/2  =  ^+i/2j+i/2Q,+i/2,7+i/2-  Finally,  the  interpolant  of  the  solution 
through  the  Gauss  points  is  a  polynomial  in  =  P^y_j  0  P^_j ; 

N-\N-\ 

Q-iXJ)  =  (37) 

1=0  7=0 

The  points  of  type  “b”  in  Fig.  7  form  the  Lobatto/Gauss  grid 
ij  =  0,1,. ..,77.  On  this  grid  are  evaluated  the  horizontal  flux  vector,  F  and 
the  metric  terms  yy  and  xy  .  The  metric  terms  are  computed  as  dy^ {X-,Y /  ^T  and 
3x^(X,,F^+i/2)/  5F.  At  points  interior  to  a  subdomain,  the  horizontal  flux  is  computed 
by 

where  Q  is  a  polynomial  of  the  type  (37)  that  passes  through  the  values 
Q(+i/2,i+i/2 !  ^i+i/2,;+i/2-  Ths  computation  of  the  flux  at  boundary  and  interface  points  is 
described  in  the  next  sub-section. 

The  vertical  flux  and  the  derivatives  yx  and  xx  are  computed  on  the 
Gauss/Lobatto  grid,  marked  by  “c”  on  Fig.  7.  The  points  on  this  grid  are 
(X;^i/22Fy),  i,j  =  0,1,..., TV  - 1.  The  metric  terms  are  computed  as  dy^ ■) !  dX  and 
dx\X,,y„Y^  )/dX.  The  vertical  flux  is  computed  at  interior  points  by 

and  at  boundary  points  as  described  in  Section  3.3. 

It  may  appear  that  to  define  quantities  on  three  different  grids  would  lead  to  a 
significantly  more  complicated  method  than  a  single  grid  Lobatto  approximation.  This 
turns  out  not  to  be  the  case.  The  definitions  of  the  fluxes  on  the  staggered  grid  by  (38) 
and  (39)  mean  that  the  reconstruction  procedure,  i.e.,  the  interpolation  of  the  solution 
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needed  to  compute  the  fluxes  at  the  Lobatto  points,  is  not  a  two-dimensional  operation. 
Rather,  it  is  the  less  expensive  sequence  of  one-dimensional  interpolations,  given  by  (11). 
The  values  of  the  solution  vector  required  to  compute  the  flux  vectors  are  actually 


and 


_  _  N-lN-l _  _ 

—  Qi+i/?..  ;+i/?A'+i/?,(-^i)^;+i/2(^>+i/2) 

i=0  j=0 

N-l _ 

“  Qi+1/2.  r+l/2^f+1/2  (■^r) 

i=0 


(40a) 


N-lN-l 


^  ^  Qi+l/2,;+l/2^i‘+l/2  (-^i+1/2  )^;+l/2  (^;) 
i=0  j=0 


since,  by  construction. 


N-l _ 

■  Qi+l/2,y+l/2^/+l/2  (Fp 

2=0 


j+\l2^  ~  ^m.j 
+1/2  (•^1+1/2)  ~  ^n.i 


(40b) 


3.3  Interface  and  boundary  treatment. 

To  describe  how  we  compute  the  interface  and  boundary  conditions  using  the 
staggered  grid  approximation,  we  will  refer  to  Fig.  8,  which  schematically  represents  four 
subdomains  and  the  locations  at  which  solution  and  flux  values  are  computed.  Only  the 
collocation  points  near  the  boundaries  are  marked.  The  circles  represent  the  solution 

values,  which  are  located  on  the  Gauss/Gauss  grid.  The  locations  of  the  horizontal  flux 
values,  F,  ,  are  represented  by  solid  squares.  The  locations  of  the  vertical  flux 

values,  marked  by  hollow  squares.  From  the  diagram,  we  see  that  along  the 

interfaces  between  subdomains  1  and  2  and  between  subdomains  3  and  4,  only  the 
horizontal  fluxes  need  to  be  computed.  Along  horizontal  interfaces,  like  those  between 
subdomains  1  and  3,  only  the  vertical  flux  needs  to  be  computed.  Because  the  grid  is 
fully  staggered,  the  coupling  is  through  subdomain  faces  only,  not  through  the  corners. 
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Fig.  8.  Diagram  of  four  subdomains  showing  locations  near  interfaces  where  solutions 
and  fluxes  are  computed.  Symbols:  %solutioi'i;  ■,F;n,G 

Fig.  8  indicates  a  significant  advantage  of  the  fully  staggered  grid  over  an 
unstaggered  grid.  In  the  unstaggered  approximation,  for  example  as  described  in  [26], 

I 

I 

special  comer  algorithms  must  be  devised  to  ensure  correct  propagation  of  waves  through 

i 

the  corners.  Each  special  case  must  be  coded  separately.  Also,  the  choice  of  bi¬ 
characteristics  that  determines  the  domains  of  dependence  becomes  more  complex  as  the 
number  of  subdomains/boundaries  that  come  together  at  a  point  increases,  making  the 

! 

derivation  of  these  special  cases  more  difficult.  The  staggered  approximation  does  not 
include  subdomain  corners,  so  conditions  do  not  have  to  be  specified  at  corner  points. 

i 

Any  number  of  subdomains  can  come  together  at  a  point  without  the  need  for  special 
point  approximations.  No  special  code  is  required  even  for  very  complex  subdomain 
topologies. 

The  interpolation  of  the  solution  by  (40)  produces  two  solution  values  at  an 
interface  point,  one  from  each  of  the  two  contributing  subdomains.  As  in  the  one 
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dimensional  case,  we  do  not  expect  these  two  values  to  coincide,  except  in  the  limit  of 
infinite  resolution.  A  single  flux  is  calculated,  as  described  for  the  one-dimensional 
problem,  except  that  we  only  consider  waves  propagating  normal  to  the  interface.  This 
normal  wave  approximation  is  common  for  finite  difference  approximations  [20]  and  has 
been  used  in  [18],  [16]  and  [4]  for  spectral  approximations.  We  note,  however,  that  other 
two-dimensional  wave  decompositions  are  possible,  like  those  surveyed  in  [32]. 

Physical  boundaries  can  be  viewed  as  interfaces  between  the  external  flow  and  the 
computational  region.  Wall  boundaries  can  be  computed  by  imposing  an  opposing  flow 
that  enforces  zero  normal  momentum  flux  across  the  interface.  Subsonic  inflow  and 
outflow  boundaries  can  be  computed  by  replacing  the  solution  that  would  have  come 
from  a  neighboring  subdomain  by  the  free-stream  values,  if  they  are  known.  If  the  full 
state  of  the  exterior  flow  is  not  known,  the  known  quantities  can  be  specified  and  the 
remaining  quantities  can  be  computed  by  a  characteristic  method.  Once  all  solution 
quantities  are  known  on  the  boundary,  the  flux  can  be  computed.  An  example  of  this 
approach  is  provided  in  Section  4.3.  Supersonic  outflow  boundaries  require  no  extra 
conditions. 

3.4  Discretization  of  the  equations. 

Once  the  fluxes  are  computed,  the  spatial  discretization  can  be  made.  From  the 
discrete  flux  values  are  defined  the  polynomials 

F(XT)  =  IE 

7Sn 

GiX,Y)  = 

i=0  j=0 

Derivatives  of  the  interpolating  polynomials  are  then  evaluated  at  the  Gauss/Gauss  grid 
points.  Like  the  reconstruction  procedure,  the  differentiation  of  (41)  can  also  be  done  as 
a  sequence  of  one-dimensional  operations 
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l/+l/2,/+l/2 


(42) 


dX 


n=0 


dG 

dr 


li+l/2,;+l/2 


m=0 


Because  both  interpolation  and  differentiation  operations  must  be  performed  at 
each  step,  the  total  work  of  the  staggered  grid  method  is  twice  that  of  a  method  that  only 
uses  the  Lobatto  grid.  The  new  method  requires  the  same  amount  of  work,  however,  as 
the  cell  averaged  method  [16].  The  reconstruction  procedure  in  two  space  dimensions  for 
the  cell  averaged  method  is  more  complex  than  in  one,  and  requires  the  same  amount  of 
work  as  both  the  interpolation  and  differentiation  operations  here. 

Finally,  from  the  definitions  (37)-(42),  the  semi-discrete  approximation  for  the 
solution  unknowns  can  be  written  as 


dQ 

_i_ 

'at 

dG 

=  0, 

(+l/2,;+l/2 

/  =  0,1,.. 

.,iV-l 

dt 

i+l/2,y+l/2 

^  dr 

;  =  0,1,.. 

.,iV-l 

Eq.  43  can  be  integrated  in  time  as  described  in  Section  2.2. 


3.S  Properties  of  the  staggered  grid  approximation. 

The  staggered  grid  approximation  is  both  conservative  and  free-stream  preserving. 
A  net  gain  or  loss  of  Q  is  determined  only  by  the  flux  through  the  exterior  boundaries.  If 
the  solution  is  constant  in  space,  then  the  solution  must  remain  constant  in  time,  even  in 
the  presence  of  a  spatially  varying  mapping. 

We  first  show  that  the  staggered  grid  approximation  is  conservative.  It  is 
sufficient  to  consider  the  four  subdomains  shown  in  Fig.  8.  Let  the  quadrature  weights 
vv,+i/2^  n;+i/2  be  defined  so  that 

1  1  N-lN-\ 

0  0  <=o  i=o 
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By  exactness  of  the  quadrature,  the  sum  of  eq.  (43)  times  Wi+in''lj+i/2  over  all  the  points 
within  a  subdomain  is 


\\^dXiY=Y,^ 

JJ  dt 


0  0 


.;=o 


N-\ 


dt 


W; 


i+1/2  Hy+lfi 


!+1/2j-H/2 


=  x 

1  1 

=  11 


0  0 


dx'^  dY 


ay 


^i+l/2^j+l/2 


(45) 


li+l/2,;+l/2 

toy 


Thus,  for  each  subdomain. 


,11  1  1 

— j  J  QdXdY  =  -  j  F(i,  y)^iy + j  F(o,  y)  jy 
^^00  0  0 

1  1 


(46) 


-|g(X,1)^X  + jG(X,0)^iX 

0  0 

When  (46)  is  summed  over  all  subdomains,  the  interior  integrals  cancel  so  that  only  the 
boundary  contributions  remain; 

d  ^  ^ 


— Z  i  1  Q'dXdY  =  j  (f'  (0,  y )  +  F'(0,  Y))dY  - 1  (F'  (1,  F)  +  F"  (1,  y))^;fy 

dt  i.=i  0  0  ^ 


(47) 


j (g' (X,0)  +  G"  (X,0)) JX  -  J  (g‘  (X,  1)  +  G'  (X,  d)  JX 

/s  A 


The  staggered  grid  approximation  is  also  free-stream  preserving,  which  means 
that  the  isoparametric  spatial  mappings  do  not  introduce  false  source  terms.  It  is 
sufficient  to  consider  the  approximation  within  one  subdomain,  since  all  derivatives  are 
computed  locally  by  subdomain.  If  we  take  F(Q)  =  G(Q)  =  1,  then  the  approximation 
(43)  becomes 


dt 


+  1 


i+l/2j+l/2 


— (yy  y^+^x) 


=  0 


(48) 
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(49) 


i+l/2J+l/2 


SO  that. 


dt 


0 


!+1/2,_;+1/2 


7=0,1,...,7V-1 


(50) 


4.  Examples 

In  this  section,  we  use  the  staggered  grid  approximation  to  compute  three  steady 
flow  problems.  The  first  problem  is  subsonic  flow  from  a  point  source,  which  has  an 
exact,  analytic  solution.  We  use  this  solution  to  show  that  exponential  convergence  is 
obtained.  The  second  problem  is  a  subsonic  flow  over  a  circular  bump  in  a  channel. 
Although  there  is  no  exact  solution  for  this  problem,  we  show  that  the  errors  due  to 
entropy  generated  along  the  curved  wall  decay  exponentially  fast.  The  final  problem 
computes  a  transonic  flow  in  an  axisymmetric  converging-diverging  nozzle.  That 
solution  is  compared  to  experimental  data. 


4.1  Subsonic  Point  Source  Flow 

As  our  first  example,  we  consider  the  flow  of  a  steady,  irrotational  flow  exiting 
from  a  point.  This  flow  can  be  solved  exactly  by  a  hodograph  transformation  [12].  The 
streamlines  are  radial,  and  level  curves  of  the  Mach  number,  pressure  and  density  are 
circles  centered  on  the  source.  We  will  compute  this  flow  in  two  geometries.  The  first 
represents  a  flow  in  an  expanding  duct,  where  two  streamlines  are  chosen  as  walls  of  the 
duct.  The  second  geometry,  a  square  with  five  circles  cut  out  of  its  interior,  is  included  to 
show  that  the  method  can  be  used  to  compute  a  flow  in  a  complex,  multiply  connected 
geometry. 

The  first  geometry  represents  steady  flow  in  an  expanding  two-dimensional  duct 
with  straight  walls.  The  lower  wall  was  chosen  to  be  the  line  y  =  0  and  the  upper  wall  was 
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the  line  y  =  x  tan(7t/6),  for  x  between  1  and  1.5.  The  exact  solution  chosen  sets  the  Mach 
number  at  the  lower  left  corner  of  the  domain  to  be  M  =  0.6. 

We  examine  solutions  for  three  subdomain  decompositions,  each  having  four 
subdomains.  Fig.  9  shows  the  three  decompositions.  In  the  first  (Grid  I),  the  subdomain 
boundaries  are  straight  lines  so  that  the  mappings  defined  by  (34)  become  bi-linear 
transformations.  The  second  and  third  decompositions  are  included  to  study  the  effect  of 
curved  subdomains.  Both  perturb  the  Grid  I  by  a  sine  wave  of  amplitude  0.1  into 
“bulging”  (Grid  H)  and  “wedging”  (Grid  III)  decompositions,  named  so  in  [36]. 


Fig.  9.  Three  subdomain  decompositions  for  the  diverging  duct  problem. 

Wall  conditions,  applied  as  described  in  the  previous  section,  are  specified  on  the 
top  and  bottom  boundaries.  The  left  boundary  is  a  subsonic  inflow  boundary.  For  that, 
we  specify  the  exact  solution  as  the  incoming  condition  for  the  Riemann  solver.  The 
right  boundary  is  a  subsonic  outflow  boundary,  and  again  the  exact  solution  is  used  to 
specify  the  external  flow.  A  perturbation  of  the  exact  solution  was  used  as  the  initial 
condition. 

Fig.  10  shows  the  computed,  steady  Mach  number  contours  for  Grid  I.  In  that 
figure,  and  in  those  following,  contour  lines  are  plotted  using  solution  values  interpolated 
from  the  Gauss  points  to  the  Lobatto  points  using  (37).  The  grids  in  Fig.  9  show  those 
Lobatto  points.  The  solutions  are  represented  interior  to  each  “cell”  bounded  by  the  grid 
lines.  The  interpolation  is  done  for  display  reasons,  since  a  plot  using  the  Gauss  points 
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would  show  gaps  between  the  subdomains,  a  result  of  the  fact  that  the  solution  is  not 
defined  at  the  interfaces.  Plotting  the  interpolant  does  give  some  visual  indication,  of  the 
size  of  the  solution  jumps  at  the  interfaces. 


The  staggered  grid  approximation  is  exponentially  convergent  for  the  point  source 
flow  in  the  duct.  For  Grids  I-III,  Fig.  1 1  shows  the  weighted  L?-  error  in  the  density  as  a 
function  of  the  polynomial  order  used  in  each  subdomain.  The  most  marked  observation 
is  that  for  this  problem,  the  error  and  the  convergence  rate  are  not  sensitive  to  the 
presence  of  curved  interfaces.  In  fact,  the  convergence  rate  for  the  “bulged” 
decomposition  is  slightly  higher  than  for  the  straight  sided  subdomains.  This  contrasts 
strongly  with  the  observations  of  [36],  which  considered  the  approximation  of  second 
order  problems.  There,  the  presence  of  even  slightly  curved  interfaces  increased  the  error 
by  orders  of  magnitude. 
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Fig.  11.  Convergence  of  the  error  for  the  three  grids  of  Fig.  9. 

As  our  last  example  of  the  point  source  flow,  we  use  the  grid  shown  in  Fig.  12. 
The  geometry,  a  square  with  five  circles  cut  out  of  its  interior,  was  chosen  to  show  that 
the  method  can  be  used  to  compute  on  a  complex,  multiply  connected  region.  Twenty- 
four  subdomains  were  used  to  cover  the  computational  domain.  Up  to  seven  subdomains 
share  a  common  corner  point  without  difficulty,  because  such  points  are  not  included  in 
the  discrete  approximation. 
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Fig.  12.  Grid  for  point  source  problem. 


The  boundary  conditions  were  chosen  so  that  the  exact  steady  solution  was  radial 
flow  with  the  point  source  at  the  center  of  the  middle  circle  of  Fig.  12.  The  center  cutout 
circle  was  specified  as  an  inflow  boundary,  with  the  conditions  chosen  so  that  the  Mach 
number  of  the  incoming  flow  was  M  =  0.6.  The  boundary  conditions  along  the  remaining 
cutout  circles  were  either  inflow  or  outflow,  depending  on  the  direction  of  the  normal 
velocity.  The  square  outer  boundary  was  an  outflow  boundary.  For  all  inflow/outflow 
boundaries,  the  exact  solution  was  used  to  provide  the  external  flow  values  required  by 
the  Riemann  solver. 

In  Fig.  13,  we  plot  the  exact  and  computed  Mach  number  contours  for  the  solution 
of  the  point  source  flow.  For  the  grid  shown  in  Fig.  12,  the  contour  lines  of  the  exact 
solution,  which  are  plotted  with  dashed  lines,  are  coincident  with  the  solid  contour  lines 
of  the  computed  solution. 
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Fig.  13.  Solution  of  the  point  source  flow  for  the  goemetry  shown  in  Fig.  12.  The  exact 
solution  is  plotted  with  dashed  lines,  the  computed  with  solid  lines. 

The  approximation  on  the  grid  of  Fig.  12  converges  exponentially.  Fig.  14  shows 
the  maximum  error  in  the  density  as  a  function  of  the  polynomial  order  in  each 
subdomain.  We  see  that  doubling  the  number  of  points  per  subdomain  causes  the  error  to 
decay  by  approximately  two  orders  of  magnitude. 

4.2  Subsonic  Flow  Over  a  Circular  Bump  in  a  Channel 

The  second  example  is  that  of  a  Mach  0.3  subsonic  flow  over  a  circular  bump  in  a 
channel.  The  geometry  and  grid  with  =  9  is  shown  in  Fig.  16.  Wall  boundaries  were 
specified  at  the  top  and  the  bottom.  At  the  left  and  right  boundaries,  the  uniform  flow 
free  stream  solution  was  specified  as  input  to  the  Riemann  solver.  Initially,  the  free 
stream  solution  was  specified  everywhere,  and  then  the  boundary  conditions  were 
imposed.  This  problem  does  not  have  an  exact  analytic  solution.  However,  since  the 
incoming  flow  was  chosen  to  be  irrotational  and  isentropic,  the  entropy  should  be  zero 
everywhere.  The  fact  that  this  is  not  the  case  can  be  the  result  of  the  spatial 
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approximation  and  to  the  normal  wave  model  used  at  the  interfaces  and  boundaries  for 
calculating  the  flux  [32], 


N 


Fig.  14.  Convergence  of  the  density  error  for  the  solution  shown  in  Fig.  13. 

Solution  contours  of  the  Mach  number  for  the  grid  shown  in  Fig.  15  are  presented 
in  Fig.  16.  The  wall  pressure  along  the  bottom,  plotted  as  the  pressure  coefficient 
Cp  =  ip-l)/  J,  is  shown  in  Fig.  17.  Finally,  a  convergence  study  of  the  entropy  errors 

as  a  function  of  N  is  shown  in  Fig.  18.  In  that  figure,  we  plot  the  maximum  value  of  the 
quantityX  =  p/ p'"  - 1,  which  should  be  zero  everywhere.  We  see  that  the  error  due  to 
entropy  generation  converges  exponentially  fast. 
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Fig.  15.  Geometry  and  grid  for  N  =  9  for  the  flow  over  a  circular  bump  in  a  channel. 


Fig.  18.  Convergence  of  the  entropy  generated  by  the  staggered  grid  approximation. 

4.3  Transonic  Flow  in  a  Converging-Diverging  Nozzle 

As  an  example  of  a  transonic  problem,  we  compute  the  flow  in  an  axisymmetric 
converging- diverging  nozzle.  We  have  chosen  the  nozzle  used  in  the  experimental 
investigation  of  Cuffel  et  al.  [13],  which  was  designed  to  show  significant  two 
dimensional  effects.  The  nozzle  consists  of  a  converging  section  with  half  angle  of  45® 
and  a  diverging  section  with  half  angle  of  15®.  The  experimental  tests  were  done  in  air 
with  a  stagnation  temperature  of  540  R  and  stagnation  pressure  of  70  psia.  The  nozzle 
geometry  and  the  grid  that  were  used  in  our  computations  are  shown  in  Fig.  19.  Note  that 
we  have  varied  the  number  of  points  per  subdomain  in  this  problem. 

To  match  the  experimental  conditions,  we  scaled  the  equations  (31)  and  (32)  by 
p  =  p*/p„,,/7  =  p*/p,„,,  where  the  represents  the  dimensional  quantity.  Under  this 
scaling,  the  temperature  and  entropy  are  T  =  T*  /  =  0.  The  initial  condition  for 

the  computation  was  the  exact  solution  of  the  quasi-one-dimensional  nozzle  that  has  same 
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area  as  the  two-dimensional  nozzle.  For  the  inflow  condition  at  the  left  boundary,  we 
specified  that  the  tangential  velocity  be  zero,  the  entropy  be  zero,  and  the  temperature  be 
unity.  At  the  right  boundary,  the  outflow  is  supersonic,  so  that  no  boundary  condition  is 
necessary. 


Since  not  all  of  the  external  flow  values  are  known  at  the  left  boundary, 
particularly  the  inflow  velocity,  it  is  not  convenient  to  use  (26)  to  impose  the  boundary 
condition.  Instead,  we  use  the  following  characteristic-like  method  that  allows  us  to 
specify  only  the  parameters  that  are  known.  The  fact  that  the  inflow  condition  sets  v  =  0 
means  that  the  flow  is  essentially  one  dimensional.  Then  we  can  write  a  left-going 
Riemann  invariant  for  the  flow.  In  terms  of  the  Mach  number,  M,  and  the  sound  speed,  a, 
that  invariant  must  satisfy 


a 


M-- 


7-1 


R~  =  n 

\omputed  ^comp. 


M. 


comp. 


7-1 


(51) 


where  the  computed  quantities  represent  values  interpolated  to  the  boundary  by  the 
reconstruction  procedure.  The  boundary  conditions  fix  the  total  temperature  and  the 
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entropy.  With  the  scaling  given  above,  this  fixes  the  total  sound  speed  at  =  Vf-  Then 
the  scaled  sound  speed  and  the  Mach  number  are  related  by 


a  - 


1  + 


(52) 


Combining  (51)  and  (52),  an  equation  for  the  Mach  number  at  the  boundary  can  be 
written 


1  + 


7-1 


M-- 


M-  ^ 


7-1 


=  R. 


'computed  * 


(53) 


Equation  (53)  can  be  written  as  a  quadratic  equation  in  the  Mach  number  and 
solved  directly.  Once  the  inflow  Mach  number  is  known,  the  sound  speed  can  be 
computed  using  (52).  From  the  Mach  number,  the  sound  speed,  tangential  velocity  and 
the  entropy,  all  remaining  variables  can  be  computed.  From  the  full  state  on  the 
boundary,  the  boundary  flux  can  be  evaluated. 

Results  computed  for  the  nozzle  are  shown  in  Figs.  20  -  23.  Contours  for  the 
pressure  are  shown  in  Fig.  20.  A  comparison  of  the  Mach  contours  and  measured  Mach 
number  in  the  neighborhood  of  the  nozzle  throat  is  shown  in  Fig.  21.  We  see  good 
agreement  between  the  computed  Mach  contours  and  the  measured  values  for  Mach 
numbers  up  to  about  1.6.  We  note  that  the  discrepancies  between  the  computed  and 
measured  Mach  numbers  are  consistent  with  the  discrepancies  observed  with  the 
solutions  of  the  inviscid  flow  solvers  reported  in  [13].  Finally,  in  Figs.  22  and  23,  we 
show  a  comparison  between  the  computed  and  measured  values  of  the  pressure  and  Mach 
numbers  along  the  upper  wall  of  the  nozzle. 
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Figure  21.  Comparison  of  computed  and  measured  Mach  contours 


Mach  Number 


Fig.  22.  Comparison  of  computed  and  measure  wall  pressure  as  a  function 
of  distance  fivm  the  nozzle  throat. 


Fig.  23.  Comparison  of  computed  and  measured  Mach  numbers  as  a 
function  of  distance  from  the  nozzle  throat. 
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5  Concluding  Remarks. 

We  have  presented  a  new,  staggered-grid  Chebyshev  spectral  multidomain 
method  for  the  solution  of  inviscid  compressible  flow  problems.  The  solutions  are 
defined  at  the  nodes  of  a  Gauss  quadrature  rule,  while  the  fluxes  are  evaluated  at  the 
nodes  of  a  Gauss-Lobatto  rule.  We  have  applied  the  method  to  one  and  two  dimensional 
problems,  but  it  should  extend  directly  to  three  dimensions. 

The  staggered  grid  multidomain  method  for  compressible  flow  problems  has 
many  desirable  features.  These  features  include 

•  Conservation. 

Mass,  momentum  and  energy  are  conserved  globally 

•  Free-stream  preservation. 

A  uniform,  steady  flow  stays  uniform  and  steady  even  for  complex 

subdomain  shapes. 

•  Temporal  accuracy. 

•  Geometric  Flexibility. 

The  domain  need  only  be  decomposable  into  quadrilaterals. 

•  Programming  simplicity. 

Corners  of  subdomains  are  not  included  as  part  of  the  approximation. 

Special  cases  do  not  need  to  be  coded.  Subdomains  have  at  most  four 

neighbors  in  two  space  dimensions  and  six  in  three  space  dimensions. 

Boundary  conditions  require  no  special  comer  treatments. 

While  the  new  method  is  flexible,  there  remain  a  few  limitations.  We  have 
assumed  that  the  approximation  is  conforming.  Thus,  subdomains  must  meet  at  a  point 
or  along  a  full  side.  Along  a  side,  the  points  at  which  the  fluxes  are  computed  must 
coincide  with  their  neighbors.  Also,  the  domain  must  be  decomposable  into 
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quadrilaterals.  We  expect  that  the  first  restriction  can  be  eliminated  by  considering  non- 
conforming  approximations  [31], 

We  have  not  considered  the  approximation  of  shocks  in  this  paper.  Methods  for 
shock  capturing  with  spectral  methods  that  have  been  proposed  to  date  include  direct 
filtering  [22],  removal  of  the  discontinuity  plus  filtering  [6],  Flux  Corrected  Transport 
methods  [15],  [16]  and  weak  filtering  with  post-processing,  [14],  [17].  Shock  fitting 
should  also  be  possible  [23],  [28].  A  thorough  study  of  the  best  options  is  beyond  the 
scope  of  this  paper. 

We  have  emphasized  steady-state  computations  in  this  paper.  However,  the 
method  is  applicable  to  time  dependent  problems.  Some  one  dimensional  examples  are 
included  here  and  in  [27]. 

The  new  method  is  a  factor  of  two  more  expensive  than  methods  that  compute 
both  the  solution  and  the  fluxes  on  the  Lobatto  grid.  This  is  due  to  the  extra  interpolation 
from  the  Gauss  points  to  the  Gauss-Lobatto  points.  FFT  techniques  cannot  be  used  to 
compute  either  the  interpolation  or  differentiation  operations  with  the  new  method.  Thus, 
the  approximation  order  within  the  subdomains  must  be  kept  low  enough  to  where  matrix 
multiplication  is  more  efficient  than  FFTs,  a  value  that  varies  from  machine  to  machine. 
However,  while  the  Lobatto  grid  methods  require  point  operations  at  subdomain  corners, 
the  new  method  requires  vector  operations  only.  We  believe  that  the  flexibility  and 
programming  ease  of  the  new  method  compensates  for  the  additional  work. 

The  new  method  differs  from  the  cell  average  method  proposed  in  [35],  [15],  and 
[16].  The  cell  average  method  requires  a  pointwise  procedure  to  be  performed  at 
subdomain  corners  that  is  not  required  by  this  method.  The  reconstruction  and 
differentiation  operations  performed  are  different.  In  one  space  dimension  the  work 

required  by  the  staggered  grid  method  is  twice  that  of  the  cell  average  method.  However, 
the  reconstruction  procedure  for  the  cell  average  method  in  two  space  dimensions 
requires  two  matrix  multiplication  operations  per  line  in  each  direction,  as  does  the 
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staggered  grid  method,  so  the  work  is  equivalent  in  two  space  dimensions.  In  three  space 
dimensions,  the  staggered  grid  method  will  require  only  two  thirds  the  work  of  the  cell 
average  method,  making  it  more  efficient  in  that  case. 

Acknowledgments 

The  first  author  would  like  to  thank  Prof.  P.L.  Roe,  Prof.  C.  -W.  Shu,  Prof.  D.  Gottlieb 
and  Dr.  D.  Sidilkover  for  many  helpful  discussions.  This  research  was  supported  in  part 
by  the  U.S.  Department  of  Energy  through  Contract  #  DE-FC05-85ER250000,  by  the 
National  Science  Foundation  through  Grant  DMS-9404322  and  by  NAS  1-19480  while 
the  first  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Science  and 
Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001. 


43 


References 


[1]  Andreassen,  0.,  and  I.  Lie, ,  “Simulation  of  Acoustical  and  Elastic  Waves  and  their 

Interaction”,  J.  Acoust.  Soc.  Am.,  95, 1994,  pp.  171-186. 

[2]  Bemardi,  C.  and  Y.  Maday,  “A  collocation  method  over  staggered  grids  for  the  stokes 

problem”,  Int.  J.  Num  Meth.  Fluids,  8,  1988,  pp.  537-557. 

[3]  Bjprhus,  M.,  “The  ODE  Eormulation  of  Hyperbolic  PDES  Discretized  by  the  Spectral 

Collocation  Method”,  SIAM  J.  Sci.  Comp.,  In  Press. 

[4]  Cai,  W.,  “High  Order  Hybrid  Numerical  Simulations  of  Two  Dimensional  Detonation 

Waves”,  AIAA  Journal,  In  Press. 

[5]  Cai,  W.,  D.  Gottlieb  and  A.  Harten,  “Cell  Averaging  Chebyshev  Methods  for 

Hyperbolic  Problems”,  Computers  and  Mathematics  with  Applications,  24,  1992. 

[6]  Cai,  W.,  and  C.W.  Shu,  “Uniform  High-Order  Spectral  Methods  for  One-  and  Two- 

Dimensional  Euler  Equations”,  J.  Comp.  Phys.,  104,  1993,  pp.  427-443. 

[7]  Canuto,C.,  M.Y.  Hussaini,  A.  Quarteroni  and  T.A.  Zang,  Spectral  Methods  in  Fluid 

Mechanics.  New  York:  Springer-Verlag,  1987. 

[8]  Carpenter,  M.H.  and  C.A.  Kennedy,  “Fourth-Order  2N-Slorage  Runge-Kutta 

Schemes”,  NASAT.M.  109112,  1994. 

[9]  Chen,  M-Q  and  C.  Chiu  ,  “Optimal  m-Stage  Runge-Kutta  Schemes  for  Steady-State 

Solutions  of  Hyperbolic  Equations”  Num..  Meth.  for  P.D.E’s,  9,  1993,  pp.  643- 

666. 

[10]  Chiu,  C., “Optimal  one-stage  and  two-stage  schemes  for  steady  solutions  of 

hyperbolic  problems”,  App/.  Num.  Math.,  11,  1993,  pp.  475-496. 

[11]  Chiu,  C.  and  D.A.  Kopriva,  “An  optimal  Runge-Kutta  Method  for  Steady-state 

Solutions  of  Hyperbolic  Systems”,  SIAMJ.  Num.  Anal.,  29,  1992,  pp.  425-438. 

[12]  Courant,  R.  and  K.O.  Friedrichs,  Supersonic  Flow  and  Shock  Waves  (Springer- 

Verlag,  New  York,  1976. 

[13]  Cuffel,  R.F.,  L.F.  Back  and  P.  F.  Massier,  “Transonic  Flowfield  in  a  Supersonic 

Nozzle  with  Small  Throat  Radius  of  Curvature”,  AIAA  J.,  7,  1969,  pp.  1364-1366. 

[14]  Don,  W.S.,  Numerical  Study  of  Pseudospectral  Methods  in  Shock  Wave 

Applications”,  J.  Comp.  Phys.,  110,  1994,  pp.  103-111. 

[15]  Giannakourous,  J.  and  G.E.  Karniadakis,  “Spectral  Element-FCT  Method  for  Scalar 

Hyperbolic  Conservation  Laws”,  Int.  J.  Num.  Meth  Fluids,  14,  1992,  p.  707. 

[16]  Giannakouros,  J.  and  G.  E.  Karniadakis  ,  “A  spectral  element-FCT  Method  for  the 

Compressible  Euler  Equations”,  J.  Comp.  Phys,  115,  1994,  pp.  65-88. 


44 


[17]  Gottlieb,  D.  and  C.-W.  Shu,  “On  the  Gibbs  Phenomenon  V:  Recovering  Exponential 

Accuracy  from  Collocation  Point  Values  of  a  Piecewise  Analytic  Function” 
ICASE  Report  94-61. 

[18]  Hanley,  P.  (1993)  “A  strategy  for  the  efficient  simulation  of  viscous  compressible 

flows  using  a  Multi-domain  pseudospectral  method”,  J.  Comp.  Phys,  108,  pp. 
153-158. 

[19]  Hesthaven,  J.S.,  “A  Stable  Penalty  Method  for  the  Compressible  Navier-Stikes 

Equations.  II.  One  Dimensional  Domain  Decomposition  Schemes.”  Submitted  to 
SIAM  J.  Sd.  Comp. 

[20]  Hirsch,  C.  (1990)  Numerical  Computation  of  Internal  and  External  Flows  Vol  II. 

John  Wiley  and  Sons,  Chichester,  England. 

[21]  Hu,  F.Q,  M.Y.  Hussaini,  and  J.  Manthey,  “Low-Dissipation  and  -Dispersion  Runge- 

Kutta  Schemes  for  Computational  Acoustics”,  ICASE  Report  94-102. 

[22]  Hussaini,  M.Y.,  D.A.  Kopriva,  M.D.  Salas  and  T.A.  Zang  "Spectral  Methods  for  the 

Euler  Equations:Part  II  -  Chebyshev  Methods  and  Shock- fitting"  AIAA  J.,  23, 
1985,  p.  234. 

[23]  Hussaini,  M.Y.,  Kopriva,  D.A.,  Salas,  M.D.  ,Zang,  T.A.  "Spectral  Methods  for  the 

Euler  Equations:  Part  I-  Fourier  Methods  and  Shock  Capturing"  AIAA  J.  ,  23, 
1985,  pp.  64-70. 

[24]  Kopriva,  D.  A.,  "A  Spectral  Multidomain  Method  for  the  Solution  of  Hyperbolic 

Systems"  Appl.  Num.  Math.  ,  2,  1986,  pp.  221-241. 

[25]  Kopriva,  D.A.,  "Solution  of  Hyperbolic  Equations  on  Complicated  Domains  with 

Patched  and  Overset  Chebyshev  Grids"  SIAM  J.  Sd.  Stat.  Comp.,  10,1989,  pp. 
120-132. 

[26]  Kopriva,  D.A.,  "Multidomain  Spectral  Solution  of  the  Euler  Gas-Dynamics 

Equations"  J.  Comp.  Phys.  ,  96 , 1991,  pp  428-450. 

[27]  Kopriva,  D.A.  and  J.H.  Kolias,  “Solution  of  Acoustic  Workshop  Problems  by  a 

Spectral  Multidomain  Method”,  to  appear  in  the  Proceedings  of  the  ICASE/LaRC 
Aeroacoustics  Workshop.  Oct.  1994. 

[28]  Kopriva,  D.A.,  T.A.  Zang,  and  M.Y.  Hussaini,  “Spectral  Methods  for  the  Euler 

Equations:  The  Blunt  Body  Problem  Revisited”,  A/AA  Journal,  9,  1991,  pp.  1458- 
1462. 

[29]  Lie,  L,  “Multidomain  Solution  of  Advection  Problems  by  Chebyshev  Spectral 

Collocation”,  J.  Sd.  Comp.,  9,  1994,  pp.  39-64. 

[30]  Liou,  M.-S.  and  C.  J.  Steffen  Jr.,  J.  Comp.  Phys,  107,  1993,  pp.  23-39. 


45 


[31]  Mavriplis,  C.A.,  “Nonconforming  Discretizations  and  a  Posteriori  Error  Estimators 

for  Adaptive  Spectral  Element  Techniques”,  Ph.D.  Dissertation,  MIT,  Feb.  1989. 

[32]  Paillere,  H.,  H.  Deconinck,  R.  Struijs,  P.  Roe,  L.  Mesaros,  and  J-D  Muller, 

“Computations  of  Inviscid  Compressible  Flows  Using  Fluctuation-Splitting  on 
Triangular  Meshes”,  AIAA  Paper  93-3301. 

[33]  Quarteroni,  A.,  "Domain  Decomposition  Methods  for  Systems  of  Conservation 

Laws:  Spectral  Collocation  Approximations",  SIAM  J.  Sci.  Statist.  Comp.,  11, 
1990,  pp.  1029-1052. 

[34]  Roe,  P.L.,  “Approximate  Riemann  Solvers,  Parameter  Vectors,  and  Difference 

Schemes”,  J.  Comp.  Phys.,  43,  1981,  pp.  357-372. 

[35]  Sidilkover,  D.  and  G.E.  Karniadakis,  “Non-Oscillatory  Spectral  Element  Chebyshev 

Method  for  Shock  Wave  Calculations”,  J.  Comp.  Phys.,  107,  1993,  pp.  10-22. 

[36]  Schneidesch,  C.  R.  and  M.O.  Deville,  “Chebyshev  Collocation  Method  and  Multi- 

domain  Decomposition  for  Navier-Stokes  Equations  in  Complex  Curved 
Geometries”,  J.  Comp.  Phys.,  106,  1993,  pp.  234-257. 

[37]  Tam,  C.K.W.,  “Overview  of  Acoustic  Workshop  Solutions”,  to  appear  in  the 

Proceedings  of  the  ICASE/LaRC  Aeroacoustics  Workshop.  Oct.  1994. 

[38]  Van  Leer,  B.,  “Flux  vector  splitting  for  the  Euler  equations”,  in.  Proc.  8th 

International  Conference  on  Numerical  Methods  in  Fluid  Dynamics, 
Berlin:  Springer  Verlag,  1982. 

[39]  Williamson,  J.H. ,  “Low  Storage  Runge-Kutta  Schemes”.  J.  Comp.  Phys.,  35,  1980, 

pp.  48-56. 


46 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reportine  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
latheriniTand  maintaining  the  data  needed,  and  completing  and  reviewingthe  collection  of  information.  Send  comments  regardingthis  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducingthis  burden,  to  Washington  Headquarters  Services,  Directorate  for  InfornrationOperatiO^  and  Reports  1215  ^ 

Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20603. 


1.  AGENCY  USE  ONLY('/.eave  blank)  2.  REPORT  DATE 

March  1995 


3.  REPORT  TYPE  AND  DATES  COVERED 

Contractor  Report 


TITLE  AND  SUBTITLE 

A  CONSERVATIVE  STAGGERED-GRID  CHEBYSHEV 
MULTIDOMAIN  METHOD  FOR  COMPRESSIBLE  FLOWS 


5.  FUNDING  NUMBERS 

C  NASl-19480 
WU  505-90-52-01 


6.  AUTHOR(S) 

David  A.  Kopriva 
John  H.  Kohas 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23681-0001 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  95-18 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

NASA  CR-195060 
ICASE  Report  No.  95-18 


11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 
Final  Report 

Submitted  to  Journal  of  Computational  Physics 


12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

U  nclassified-U  nlimited 


12b.  DISTRIBUTION  CODE 


Subject  Category  64 


13.  ABSTRACT  (Maximum  200  words) 

We  present  a  new  multidomain  spectral  collocation  method  that  uses  staggered  grids  for  the  solution  of  compressible 
flow  problems.  The  solution  unknowns  are  defined  at  the  nodes  of  a  Gauss  (quadrature  rule.  The  fluxes  are  evaluated 
at  the  nodes  of  a  Gauss-Lobatto  rule.  The  method  is  conservative,  free-stream  preserving  and  exponentially  accurate. 
A  significant  advantage  of  the  method  is  that  subdomain  corners  are  not  included  in  the  approximation,  making 
solutions  in  complex  geometries  easier  to  compute. 


14.  SUBJECT  TERMS  15.  NUMBER  OF  PAGES 

Spectral  Methods;  Domain  Decomposition;  Gas  Dynamics 

16.  PRICE  CODE 

_ A03 _ 


18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION 
OF  THIS  PAGE  OF  ABSTRACT  OF  ABSTRACT 

Unclassified 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

U  nclassified 


NSN  7540-01-280-5500 


Standard  Form  298(Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


